/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */


/*! \file FinitieDifferenceDelta.h
    \brief Multidimensional generalization of finite difference Laplacian operators.

	Multidimensional generalization of the classical Laplacian difference operator
	\f[
			\delta_{x\overline{x}} u_i=\frac{u_{i-1}+2u_i-u_{i-1}}{h^2}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/
#if !defined(_u2_Delta_)
#define _u2_Delta_
#include <src/math/Operators/Delta.h>

/*! 
    \brief \f$u^p\f$ function.

	Common behavior for \f$u^p\f$ function, i.e. for a field \f$\{u_i\}_i\f$ it defines recursively
	\f[
			\{u^p(u)\}_i=\{u^p_i\}_i
	\f]
*/
template <class A,class R, OS_int P>
struct u_p
	: public OS_Generic_Algebra<A , R> 
{
	enum{Dim = A::Dim};
	typedef u_p													type;
	typedef OS_Generic_Algebra<A, R>							Base;
	typedef typename Base::argument_base_type					argument_base_type;
	typedef Base												Base_Cast;
	typedef Null												interface_type;

	u_p() {};
	u_p(const Base &from) : Base(from)
	{};
	virtual ~u_p(){};
	virtual smart_ptr_type Create_Instance() const {
		return Create_Instance_helper<type>();
	};
	virtual result_type operator()(typename argument_type from) const 
	{
		if (first_)
			return pow( (*first_)(from),P);
		return pow(from,P);
	};
};

/*! 
    \brief Common multidimensional behavior for all finite difference Laplacian operators.

	Common multidimensional behavior for all finite difference Laplacian operators, i.e.
	\f[
			|u|^2 \delta_{x\overline{x}} u_i=|u|^2 \frac{u_{i-1}+2u_i-u_{i-1}}{h^2}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/

template <class A, std::size_t = A::Dim>
class u_2_Delta 
	: public OS_Generic_Algebra<A , A> 
	, public virtual graph_base<A>

{
public :
//! return " u^2 \Delta u"\f$.
		static std::string get_name(){return " u^2 Delta u";};
		enum:std::size_t												{Order = 2};
		typedef OTSCentered												Centered_Type;
		typedef typename A::Boundary_type								Boundary_type;
		
		typedef OS_Generic_Algebra<A , A>								Base;
		typedef  Base													Base_Cast;
		typedef u_2_Delta												type;
		typedef Finite_Diff_Operator									Family_Operator;


		typedef u_p<A,A,2>												u_2_type;
		typedef Delta<A, A::Dim, smart_ptr<A> >							Delta_p_t;
		typedef typename Delta_p_t::graph_type							graph_type;
		typedef typename Delta_p_t::weight_set_type					weight_set_type;

		//! return " \delta_{\overline{x}x}u"\f$.
		type( ) {};
		type(const type & data) : Delta_(data.Delta_), Base(data) {};
		type(const weight_set_type & data) :  Delta_(data) {};
		type(const Base_Cast &first) : Base(first)
		{};
		virtual ~u_2_Delta(){};
		type & operator = (const type & from) {
			Base::operator =(from);
			Delta_ = from.Delta_;
			return *this;
		};

		virtual const graph_type& get_graph() const {
			return Delta_.get_graph();
		}
		virtual void update() {
			Base::update();
			Delta_.set_state_variable(get_state_variable());
		};
		virtual result_type operator()(typename argument_type from) const 
		{
			return u_2(from)^Delta_(from);
		};
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
		Delta_p_t	Delta_;
		u_2_type u_2;
};
/*! \brief Bind operator for operator u_2_Delta

	This is the time Bind operator of Laplacian difference operator u2_Delta_traits, that is u2_Delta_t(t,u)=u_2_Delta(u)
*/
template <class data, typename Time = OS_double >
struct u_2_Delta_t : public BindTimeOperator<u_2_Delta<data, data::Dim> ,Time > {

	typedef u_2_Delta_t													time_operator_type;
	typedef BindTimeOperator<time_operator_type ,Time >					Base;
	u_2_Delta_t(){};
};



#endif